'''
This tests the OLS Model based on the Longley dataset and known results.

References
-----------

Longley, James W.  1967  "An Appraisal of Least Squares Programs for the 
Electronic Computer from the Point of View of the User"  Journal of the 
American Statistical Association.  Vol. 62, No. 318, 819 - 841.

http://www.stanford.edu/~clint/bench/

http://www.itl.nist.gov/div898/strd/

'''
from numpy.testing import *
import numpy as np
import os
import nipy.fixes.scipy.stats.models as SSM

import rpy
from rmodelwrap import RModel

filename=os.path.join(os.path.dirname(os.path.abspath(__file__)), "longley_data")
data=np.loadtxt(filename, dtype=np.float, skiprows=1, usecols=(1,2,3,4,5,6,7))
y=data[:,0]
x=data[:,1:]
# must add a constant, since it still depends on formula framework
x=np.hstack((np.ones((len(x),1)),x))
ntis_longley=(-3482258.63459582,15.0618722713733,-0.358191792925910E-01,-2.02022980381683, -1.03322686717359, -0.511041056535807E-01,1829.15146461355) # keep with the data?

def test_OLSModel():
    res=SSM.regression.OLSModel(x).fit(y)
    assert_almost_equal(res.beta, ntis_longley, 5)  # how much precision do we need?  
    # I'm sure precision machine and somewhat algorithm specific.
    # interesting to compare
    s_xx=np.linalg.inv(np.dot(x.T,x))
    s_xy=np.dot(x.T,y)
    assert_almost_equal(np.dot(s_xx,s_xy),ntis_longley,3)
    rlm_res = RModel(y, x, rpy.r.lm)
    assert_almost_equal(rlm_res.beta,ntis_longley,8)    # should we at least strive for R precision?
     

if __name__=="__main__":
    run_module_suite()


